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ANALYSIS OF SHELL-TYPE STRUCTURES SUBJECTED TO 
TIME-DEPENDENT MECHANICAL AND THERMAL LOADING 
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INTRODUCTION 

The objective of the present research is to develop a 
general mathematical model and solution methodologies for analyz- 
ing structural response of thin, metallic shell-type structures 
under large transient, cyclic or static thermomechanical loads. 
Among the system responses, which are associated with these load 
conditions, are thermal buckling, creep buckling and ratcheting. 
Thus, geometric as well as material-type nonlinearities (of high 
order) can be anticipated and must be considered in the develop- 
ment of the mathematical model. Furthermore, this must also be 
accommodated in the solution procedures. 


SUMMARY OF PROGRESS 

The progress to date has been elaborated upon in an interim 
scientific report submitted to the sponsor during the summer of 
1986, and in a series of semiannual progress reports. The most 
recent of these is dated October, 1987. 


A complete true ab-inito rate theory of kinematics and kine- 
tics for continuum and curved thin structures, without any res- 
triction on the magnitude of the strains or the deformation, was 
formulated. The time dependence and large strain behavior are 
incorporated through the introduction of the time rates of the 


metric and curvature in two coordinate systems; a fixed (spatial) 
one and a convected (material) coordinate system. The relations 
between the time derivative and the covariant derivatives (gra- 
dients) have been developed for curved space and motion , so that 
the velocity components supply the connection between the equa- 
tions of motion and the time rate of change of the metric and 
curvature tenors. 

The metric tensor (time rate of change) in the convected ma- 
terial coordinate system is linearly decomposed into elastic and 
plastic parts. In this formulation, a yield function is assumed, 
which is dependent on the rate of change of stress, metric, tem- 
perature, and a set of internal variables. Moreover, a hypoelas- 
tic law was chosen to describe the thermoelastic part of the de- 
formation. 

A time and temperature dependent viscoplastic model was for- 
mulated in this convected material system to account for finite 
strains and rotations. The history and temperature dependence 
were incorporated through the introduction of internal variables. 
The choice of these variables, as well as their evolution, was 
motivated by phenomenological thermodynamic considerations. 

The nonisothermal elastic-viscoplastic deformation process 
was described completely by "thermodynamic state" equations. Most 
investigators (in the area of viscoplasticity) employ plastic 
strains as state variables. Our study shows that, in general, use 
of plastic strains as state variables may lead to inconsistencies 
with regard to thermodynamic considerations. Furthermore, the 
approach and formulation employed by all previous investigators 
lead to the condition that all plastic work is completely dis- 
sipated. This, however, is in contradiction with experimental 
evidence, from which it emerges that part of the plastic work is 
used for producing residual stresses in the lattice, which, when 
phenomenologically considered, causes hardening. Both limitations 
are not present in our formulation, because of the inclusion of 
the "thermodynamic state" equations. 

The obtained complete rate field equations consist of the 
principles of the rate of the virtual power and the rate of 
conservation of energy, of the constitutive relations, and of 
boundary and initial conditions. These formulations provide a 
sound basis for the formulation of the adopted finite element 
solution procedures. 

The derived shell theory, in the least restricted form, 
before any simplifying assumptions are imposed, has the following 
desirable features: 

(a) The two-dimensional, impulse-integral form of the equations 
of motion and the Second Law of Thermodynamics ( Claus ius- 
Duhem inequality) for a shell follow naturally and exactly 
from their three-dimensional counterparts. 

(b) Unique and concrete definitions of shell variables such as 
stress resultants and couples, rate of deformation, spin and 


entropy resultants can be obtained in terms of weighted 
integrals of the three-dimensional quantities through the 
thickness. 

(c) There are no series expansions in the thickness direction. 

(d) There is no need for making use of the Kirchhoff Hypotheses 
in the kinematics. 

(e) All approximations can be postponed until the descretization 
process of the integral forms of the First Law of Ther- 
modynamics 

(f) A by-product of the descent from three-dimensional theory 
is that the two-dimensional temperature field (that emerges) 
is not a through-the-thickness average, but a true point by 
point distribution. This is contrary to what one finds in 
the literature concerning thermal stresses in the shell. 

To develop geometrically nonlinear, doubly curved finite 
shell elements the basic equations of nonlinear shell theories 
have to be transferred into the finite element model. As these 
equations in general are written in tensor notation, their 
implementation into the finite element matrix formulation 
requires considerable effort. 

The nonlinear element matrices are directly derived from the 
incrementally formulated nonlinear shell equations, by using a 
tensor-oriented procedure. The classical thin shell theory based 
on the Kirchoff-Love hypotheses (Formulation D) was employed for 
this purpose. For this formulation, we are using the '•natural” 
degrees of freedom per mid-surface shell node: three incremental 
velocities and the rates of rotations about the material coor- 
dinates in a mixed form. 

A description of the developed element and related finite 
element code are given in Appendix A. This exposition provides 
information concerning the formulation, the finite element and 
how it is employed in the solution of shell-like configurations. 
A complete description including program flow chart, listing and 
instructions to the user is been finalized and it will be sent to 
the sponsor, separately. 

The quasi-linear nature of the principle of the rate of 
virtual power suggests the adoption of an incremental approach to 
numerical integration with respect to time. The availability of 
the field formulation provides assurance of the completeness of 
the incremental equations and allows the use of any convenient 
procedure for spatial integration over the domain V. In the 
present instance, the choice has been made in favor of a simple 
first order expansion in time for the construction of incremental 
solutions from the results of finite element spatial integration 
of the governing equations. 

The procedure employed permits the rates of the field formu- 
lation to be interpreted as increments in the numerical solution. 
This is particularly convenient for the construction of incremen- 
tal boundary condition histories. 


Even under the condition of static external loads and slowly 
growing creep effects, the presence of snap-through buckling 
makes the inertia effects significant. In dynamic analyses, the 
applied body forces include inertia forces. Assuming that the 
mass of the body considered is preserved, the mass matrix can be 
evaluated prior to the time integration using the initial 
conf igurat ion . 

Finite element solution of any boundary-value problem 
involves the solution of the equilibrium equations (global) 
together with the constitutive equations (local) . Both sets of 
equations are solved simultaneously in a step by step manner. The 
incremental form of the global and local equations can be 
achieved by taking the integration over the incremental time step 
t=tj + ^-tj. The rectangular rule has been applied to execute the 
resulting time integration. 

Clearly, the numerical solution involves iteration. A simp- 
lified version of the Riks-Wempner constant-arc-length method has 
been utilized. This iteration procedure which is a generalization 
of the displacement control method also allows to trace the non- 
linear response beyond bifurcation points. In contrast to the 
conventional Newton-Raphson techniques, the iteration of the 
method takes place in the velocity and load rate space. The load 
step of the first solution in each increment is limited by 
controlling the length ds of the tangent. Either the length is 
kept constant in each step or it is adapted to the characteris- 
tics of the solution. In each step the triangular-size stiffness 
matrix has to be checked for negative diagonal terms, indicating 
that a critical point is reached. 

One of the most challenging aspects of finite strain formu- 
lations is to locate analytical solution with which to compare 
the proposed formulation. Typically, as a first problem, a large 
strain uniaxial test case was analyzed. The case considered 
examines the rate-dependent plastic response of a bar to a 
deformation history that includes segments of loading, unloading, 
and reloading, each occurring at varying strain and temperature 
rates. Moreover, it was shown that the proposed formulation 
generates no strain energy under a pure rigid body rotation. 
These are surely important demonstration but they only represent 
a partial test because the principal stretch directions remain 
constant. Finally, a problem which was discussed by Nagtegaal and 
de Jong, and others too, as a problem which demonstrates limita- 
tions of the constitutive models in many strain formulation, is 
the Couette flow problem. This problem is solved as a third 
example. The results of these test problems show that: 

- The formulation can accommodate very large strains and 
rotations . 

- The formulation does not display the oscillatory beha- 
vior in the stresses of the Couette flow problem. 



The model incorporates the simplification associated 
with rate-insensitive elastic response without losing 
the ability to model rate temperature dependent yield 
strength and plasticity. 

The problem of buckling of shallow arches under transient 
thermomechanical load was investigated next. The analysis was 
performed with the aid of 24 paralinear isoparametric elements. 
The paralinear isoparametric element is such that the thickness 
is small compared to other dimensions. The characteristics of the 
element are defined by the geometry and interpolation functions, 
which are linear in the thickness direction and parabolic in the 
longitudinal direction. Consequently, the element allows for 
shear strain energy since normals to a mid-surface before 
deformation remain straight, but not necessarily normal to the 
midsurface after deformation. 

The developed solution scheme is capable of predicting 
response which includes pre- and post-buckling with thermal creep 
and plastic effects. The solution procedure was demonstrated 
through several examples which include both creep and snap- 
through behavior. 

The last set of problems which are under investigation 
consists of creep or thermal buckling, with plastic effects, of 
shells of revolution. 

In addition, following a more traditional approach, a method 
was developed for bounding the response (solution) of bars and 
beams of (linear) viscoelastic material behavior, based on 
nonlinear kinematic relations. 

In connection with this progress to date, two papers were 
published by the AIAA Journal in 1986. Moreover, a paper en- 
titled, "Non-Isothermal Elastoviscoplastic Analysis of Planar 
Curved Beams" was presented at the 3rd Symposium on Nonlinear 
Constitutive Relations for High-Temperature Applications, held at 
the University of Akron, on June 11-13, 1986. A descriptive 
abstract of this paper was published in the meeting proceedings 
and the full paper will appear in a special publication. Copies 
of the above have been sent to the sponsor. 

In addition, the two papers presented at the 28th AIAA/- 
ASME/ASCE/AHS SDM Conference and published in the proceedings of 
this conference, and the one paper presented at the 21st Con- 
ference of the Israel Society of Mechanical Engineers (keynote 
lecture by Dr. R. Riff) been accepted for publication; the first 
two by the AIAA Journal and the last one by the International 
Journal of Computers and Structures. These three papers deal with 
applications to snap-through and creep buckling of bars and 
arches, and to two-dimensional problems in extension and shear. 
Most of this work was also presented at the NASA-Lewis Conference 
on Structural Integrity and Durability of Reusable Space Prop- 



ulsion Systems on May 1987 in Cleveland. Copies of these papers 
have been sent to the Sponsor. 

In connection with the more traditional approach a paper 
accepted for presentation and for publication in the Proceedings 
of the special Symposium on Constitutive Equations at the ASME 
Winter Annual Meeting, Chicago, IL. , November 28 - December 2, 
1988. The title of the paper is "Creep Analysis of Beams and 
Arches Based on a Hereditary Visco- Elastic- Plastic Constitutive 
Law" and a copy of the paper is attached herewith. 

Moreover, two papers are presently in preparation, and will 
be submitted for publication during the month of May. One deals 
with applications to shell configurations, and it is titled 
"Analysis of Shell Structures Subjected to Large Mechanical and 
Thermal Loadings" . The second deals with the time- dependent 
response of curved structural elements, and it is titled "The 
Dynamic Aspects of Thermo- Elasto- Viscoplastic Snap- Through and 
Creep Phenomena" . Copies of these papers will be sent to the 
sponsor with the next semi- annual report. 

The study was initially limited only to the simplest of the 
developed shell theory formulations (Formulation D) . 

The various shell theory approximations (formulations) are 
based on the use of certain simplifying assumptions regarding the 
geometry and kinematics of the shell configuration. 

These are: 

Assumption I : The material points which are on the normal to 
the reference surface before deformation will be on the same 
normal after deformation. 

Assumption II : The shell is sufficiently thin so that we can 
assume linear dependence of the position of any material 
point (in the deformed state) to the normal (to the refe- 
rence surface) coordinate (in the deformed state) . The 
linear dependence can easily be changed to parabolic, cubic, 
or any desired degree of approximation. 

Assumption III : The rate of change of the velocity gradients 
with respect to in-plan coordinates on the two bounding 
shell surfaces is negligibly small. 

Assumption IV : The rate of change of the distance of a 

material from the reference surface is negligibly small. 

On the basis of the above four simplifying assumptions, 
several formulations result, for the analysis of thin shells. 
These formulations are denoted below by capital letters. 

Formulation A : This formulation makes use of Assumption 

I, only. 

Formulation B : This formulation employs Assumptions I and 

II. 

Formulation C : This formulation employs Assumption I, II and 



III. 

Formulation D : This is the classical thin shell theory based 

on the Kirchhof f-Love hypotheses of Assump- 
tions I, II, III, IV, as applied to large 
deformation theory. 

These formulation are arranged in such a manner that we move from 
the least restrictive (A) to the most restrictive (D) . 

In addition to this a fifth formulation (E) can easily be 
devised and this formulation in terms of order of restriction is 
similar to Formulation A. Formulation E makes use of Assumption 
II only. 


FUTURE TASKS 


The main thrust of the additional tasks is to develop a 
finite element and select a code, which will be made available to 
all users and which will be based on the most general (but 
practical) nonlinear shell formulation possible and nonlinear 
constitutive relations to predict the response of shell-like 
structures, when subjected to time-dependent thermomechanical 
loads with large excursions. 


In order to accomplish the overall objective, as stated a 
number of steps or tasks need be completed. These are: 


Step 1 : 


Step 2 : 
Step 3 : 


Step 4 : 


Step 5 : 
Step 6 : 


Formulation the principle of the rate of virtual 
power for the less restricted shell formulation 
(C, B, E, A). 

Couple the general constitutive relations to shell 
kinematics for the aforementioned formulations. 
Before developing a finite element, and before 
numerically testing and comparing the merits of 
the various formulations, some theoretical 
evaluation of the degree of approximation will be 
attempted . 

Develop a finite element for the chosen least 
restrictive shell formulation. In order to 
accomplish this some novel ideas and procedures 
will be employed. These include a tensor-oriented 
procedure for obtaining the element, and a 
discretization process that involves not only the 
nodal degrees of freedom but also an approximation 
to the ’’shift" tensors. 

Incorporate the developed finite element and 
related constitutive relations into a FEM code. 
Numerically test and evaluate the developed finite 
element procedure. 


It is estimated that this entire research effort will 
require a minimum of three years to bring to fruition. Some of 
the tasks (steps 1 and 2) have already been started. 


APPENDIX A - THE SHFIT/T. KT.EMKNT 


A brief description highlighting the essential features of 
the shell element development and the related code used in this 
work is given here. 

In order to derive discrete algorithm based on the finite 
element displacement method we approximate the velocity field by 
index-oriented notation, which allows the separate representation 
of the shape functions (the specific expression depends on the 
decided upon degress - of - freedom, Lagrangian, Hermitian, etc.) 
for the tangential velocities and for the normal velocity 0* . 

\ji= \ t £ V/\ - v* 

< v* * v* 

Upper indices imply the columns, lower indices the rows of a 
matrix expression, and the summation is carried out spanning the 
number of degrees-of- freedom. V** and represent, therefore, the 
vector of nodal velocities by the row and by the column respec- 
tively. We get the shape functions for the partial derivatives 
of the velocity shape functions v* . * : 


vr^*p Vm 

(3) 


(4) 


The main idea of this formulation is the development of 
shape functions for further mechanical and thermal variables by 
the application of well-known tensor procedure on the basic shape 
functions (3) and (4). Taking, for example, the covariant 
derivative 

(5) 

and inserting (3) and (1) , we can define the shape function of 
the covariant derivative: 

M ^ 'T T 

vT<*} f ^ ~ )j Va = VJVp. ** (6) 

—y — " 

In the same way we receive the shape function of the 
internal variables, for example the rate of deformation and the 
spin tensors: 

^ y ; 



( 1 ) 

( 2 ) 


( 7 ) 


( 8 ) 


~ ~ 4 'F>>. Up 4 + ^2> k f>) 

~ (^i)^ 4 ^5;* kj 4 k2jp- VTj kid. \>^2 ^ = 

v — 

-<A>fi f )Vi 

The same procedure is now applied to the shift tensor 


/r';*U?*$oW fj + fj J? (9) 

which is responsible for the "exact" distributions over the 
thickness. 

Finally the shape functions of the internal forces and 
temperature variables can be derived from the shape functions of 
the rate of deformation and spin tensor via the constitutive 
relations; for example; 

= v DH*^* VT*?*), Va* Y/a ( 10 ) 

* ■ ■ ' * y " — . f 

All these expression are now substituted into the rate of 
the first law of thermodynamics to obtain the element "stiffness" 
equations. The developed element matrices are implemented into 
the global relation of the complete shell structure by standard 
assemblage process considering incidence and boundary conditions. 

The present curvilinear formulation of the element enables 
the precise description of geometry, external loads and tempera- 
tures and the fulfillment of the convergence criteria, while the 
rigid body motion condition can only be satisfied in an approxi- 
mate manner. The tensor oriented formulation renders the 
optional use of various shape functions for the tangential 
velocities v and the normal velocity v 3 . 

The shell element which have been used up to today is based 
on the bicubic Hermite polynomial with 4x12 generalized velocit- 
ies and 4 temperatures. Numerical integration spanning the 
element domain was applied (16 points of integration) , whereby 
area and boundary integrals were replaced by double integration 
with respect to the curvilinear O'* - coordinates. 
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FLOW CHART FOR THE SET-UP OF THE ELEMENT MATRICES 















CREEP ANALYSIS OF BEAMS AND ARCHES BASED ON A 
HEREDITARY VISCO-ELASTIC-PLASTIC CONSTITUTIVE LAW 

by 

John M. Stubstad* and George J. Simitses** 

Georgia Institute of Technology 
Atlanta, Georgia 

ABSTRACT 

In both the nuclear power and aerospace industries, it is not uncommon to find metallic structural 
elements which are subjected to significant thermal excursions from an elevated ambient temperature. 
Depending upon the nominal ambient temperature and the magnitude of the thermal cycle, this can 
cause the time-dependent properties of the material to change radically during each duty cycle. At 
the lower temperatures of the duty cycle, the material may behave primarily in a viscoelastic manner. 
However, as the temperature or mechanical loading increases, plastic-type behavior may predominate. 
These modifications of the behavior of the material substantially complicate the analysis of the problem, 
thereby making prediction of the ultimate useful life of such an element exceedingly difficult. 

The paper presents an analytic study of planar beams and arches subjected to significant thermal 
cycling from ambient temperatures up to 600 P C. In the study, a recently developed unified nonlin- 
ear hereditary type of viscoelastoplasiic constitutive law is employed to characterize the time- and 
temperature-dependent properties of a typical aerospace alloy, Hastelloy X. 

‘The results from this work demonstrate that a strong interaction exists between the back stress 
variable of this particular constitutive law and the time-dependent stress distribution produced by the 
geometry of the deformation. Effectively, this interaction tends to control, in a highly nonliuear manner, 
the creep : ratchetting response* of the beam and the arch. An unexpected consequence of this is that 
temperature gradients in the thickness direction, a factor normally neglected in most studies, tends to 
exert an important influence on the response during thermal cycling. 

NOMENCLATURE 

a / load eccentricity 

A cross sectional area of beam or arch 

’ b width of beam or arch 

• dj - . inelastic strain tensor * 

d deformation vector for points-on the centrcidal axis 

E Young's Modulus 

J?o, E\ y and E 2 zero, first and second moments of the elastic modulus across a 

cross section, respectively 

• currently with the Aerospace Corporation; Member of ASME 

•• Professor of Aerospace Engineering; Member of ASME 


g a , G a base vectors for undeformed and deformed configurations, respectively 

9ap G a p metric components of the undeformed and deformed configurations, respectively 
h depth of beam or arch 
k K+ 

K drag stress 

K i , Ki , n t thru n 7 constitutive law constants 
m and n constitutive law exponents 
M and N moment and force resultants, respectively 
M e and N c creep strain moment and force pseudo-resultants, respectively 
M* and thermal strain moment and force pseudo-resultants, respectively 

p pressure load 

P and V axial and transverse force resultants acting on the cross section 
Q(t) and T(t) constitutive law functions 

r and R position vectors in undeformed and deformed configurations, respectively 
s and 7) coordinates along length and depth directions, respectively 
s ij components of the deviator stress tensor 

t, n and k triad of unit vectors for the undeformed configuration 

T, N and K triad of unit vectors for the deformed configuration 

i . time 

u and w axial and transverse displacement of the centroidal axis, respectively 
a coefficient of thermal expansion 

e Q centroidal axis strain 

€{j strain tensor component 

0 change in temperature 

k initial curvature of the arch 

A and fi Lame constants 

o-ij stress tensor component 

<f> angle of rotation of cross section 

O t; backsiress 

INTRODUCTION 

It is well known that' metal alloys can undergo transitions in behavior as temperature increases. 
Commonly, for loading substantially below yield, the elastic response observed at room temperature gen- 
erally gives way to a time-dependent viscoelastic response at somewhat elevated temperatures. Further 
increases in temperature, however, introduces the potential for sudden “rapid” or plastic type deforma- 
tion. Such transitions create the potential for significantly shortening the useful life of the structural 
element and ultimately generate the possibility for a sudden unanticipated failure. Consequently, for 
many years the aerospace and nuclear power industries, where elevated temperature operating environ- 


ments abound, have had a continuing interest in predicting the behavior of metallic structural elements 
subjected to such conditions. 

Early investigators, such as Hilton and others [1-3], generally focused their attention on the 
behavior of structural components subjected to conditions of constant load at constant uniform elevated 
temperatures. Many of these studies employed simplified analyses to improve mathematical tractability. 
Additionally, “experimentally based equation of state” type constitutive laws were often used to express 
the nonlinear elevated temperature time-dependent behavior of the material. A summary of many of 
methods developed and key findings obtained from these studies is available in Hoff [4]. 

These efforts provided answers to many of the questions regarding elevated temperature creep 
buckling. However, they were not able to satisfactorily describe the creep ratchetting behavior resulting 
from elevated temperature thermal cycling. Consequently, it was not until the efforts of Miller [5] and 
Edmunds and Beer [6], both of whom considered nuclear pressure vessels, that this particular form 
of behavior was specifically addressed. Subsequently, Bree [7,8] investigated the basic factors which 
determine when this type of response could occur. Of course, these early studies lead to many subsequent 
efforts which explored various aspects of the phenomena. Some of these studies relied on experimental 
investigations while others employed analytic approaches. In this regard, the works of Conway et. al. [9], 
Corum [10,11] and Mukhejeree, Kumar and Chang [12] provide a representative sampling of the types 
of studies which were conducted. 

However, the greatest concentration of effort related to this topic area has been the associated 
effort directed toward improving the capability to predict the elevated temperature behavior of metals. 
The contributions of Hart [13,14], Pointer and Leckie [15], Pugh [16-18], Krempl [19], and Walker and 
Krempl [20], to name of few, provide a dramatic illustration of the intensity of these efforts to develop 
and employ advanced constitutive models. Yet, as pointed out by Corum and Sartory [21], an equation of 
state approach to constitutive modeling is still generally used in design situations. However, as Pugh [18] 
has noted, the newer types of unified constitutive laws, where inelastic strain is not divided into distinct 
creep and plasticity components, can provide an alternate approach. 

Consequently, one of the aims of this study is to examine the use of a typical unified constitutive 
model in an analysis of the behavior of structural elements subject to thermal cycling from an elevated 
ambient temperature. The specific law which was selected is one developed by Walker [22] to model 
the time- and temperature-dependent behavior of Hastelloy X, an alloy routinely used in the aeropsace 
industry. 

The results of this study indicate that, for this particular form of constitutive law, an implicit 
interaction exists between the -actual stress in the member and* the backstress variable of the constitutive 
law. This interaction is shown to strongly influence the ultimate response of the element. This result is 
significant for two reasons. First, because saturation of the backstress of this constitutive law can lead 
to plastic response, which signifies that the ultimate* reliability of any predicted results rests strongly 
upon the accuracy with which the backstress growth law parameters have been determined. 

Additionally, this particular aspect produces the rather interesting result that, due to the strong 
temperature dependence of the material constants of the law, temperature variations in the thickness 
direction greatly influence predicted response. This is a significant for any analysis since the influences 


of such temperature variations are generally neglected in most studies. 


MATHEMATICAL FORMULATION 

To focus principally upon the interaction between the response of the structural element and 
the prediction of thermal dependence of the material, the problem is formulated within the context of 
a simplified beam theory. Consequently, it is assumed that the beam or arch deforms in accordance 
with the Euler-Bemoulli hypotheses. As such, cross sectional planes normal to the centroidal axis in 
the undeformed geometry are assumed to remain plane and normal in the deformed state. Similarly, 
extensional straining in the thickness and depth directions are neglected. Therefore, based on the 
geometry illustrated in Fig. 1, the position vectors r and ft, where 


r = r Q -J- Tjn and R = r a + d + 77 N (1) 

are employed to locate a typical point on an arbitrary cross section in the undeformed and deformed 
configurations, respectively. Note that r) represents the coordinate in the normal direction. Also, lower 
case and upper case symbols are employed to denote quantities referred to the undeformed and deformed 
configurations, respectively. 

Base vectors for the reference and current state, g a and G a , respectively, are defined by 


g a = — and G a = — where a = a, ij ' (2) 

8a da 

Consequently, the deformation vector which translates a typical point from the undeformed to deformed 
configurations, denoted as d, can be expressed as 

d 4* = (u + ??sin^)t -f (w -f r}cos<f>)n _ (3) 

where u and w represent axial and transverse displacement functions for points on the centroidal axis, 
respectively. Thus, upon substituting Eqns. (l) and (3) into (2), differentiating and employing the 
Fernet-Serret formulae and the strain definition, 

( 4 ) 

(5a) 


7a/3 = - 

yields the strain expressions 


where - g a p = g a * g /3 , G a p = G a • Gp 


7.. = ^{(1 + + Tjfecos^) 1 + - HU - T)k sin <f>) 2 - (1 + 


8b 


7*«j — 7«m = J 1 + KW + & cos 4) sin <6 + ( - *)k sin 4>) cos if>) 


and 


(55) 


7nn = ^{sinV + cos 1 ^- l} = ° (5c) 

Note that in these expressions, k = k + d4>/ds y where k denotes the initial constant curvature of the 
arch and <{> represents the angle of rotation of the cross section. 



From the Euler* Bernoulli hypotheses, the shear strains must vanish. Consequently, from this 
requirement Eqn. (5b) yields 


tan 4> — 


,dw . 
- *“> 

( i + _ + Kw) 


Therefore, employing the definition 


(6) 


„ ^ / Qu ,du) X 2 

l + 2€ 0 = (1+ + 


ds 


and Eqn. (6) it is relatively easy to show that 


sin<£ 


,dw . 

-\T7~ KU ) 


_ Ids 


and 


cos d> — 


du 

1+ — + KW 
03 


(7) 


( 8 ) 


( 9 ) 


y/l + 2e 0 y / 1 + 2e 0 

Thus, expanding Eqn. (5a) and employing Eqns. (7) and (8) yields 

( , 8<f> ^ 2 f 1 8<f> x 8(f) 

7** =€o + + 2 € 0 ~~K+ T2e 0 — } + 9 {* + 2~dj}j7 

From the above it is clear that e a represents the strain induced along the centroidal axis. Since, for a 
thin arch or beam, the last term of Eqn. (9) should be small in comparison to the others, it may be 
neglected. Similarly, additional simplifications may be obtained for the case where the centroidal axis 
strain is sufficiently small so that it may be neglected in comparison to one. Based on these assumptions^ 
Eqn. (9) simplifies to the standard form . 
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Equilibrium equations are obtained through application of the Principle of Virtual Work, Stress 
resultants, N and M, are defined such that 
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Employing the deformed configuration geometry illustrated in Fig. 2, it is noted that 
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where 6<f> and 6e denote the incremental increases in rotation and centroidal axis strain resulting from 
the deformed configuration axial and transverse displacements, V\ and v 2) respectively. 

Consequently, from the Principle of Virtual Work, 
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Expressing the first term on the right hand side in terms of an integral over the length and then combining 
that result with the work term yields 
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Therefore, employing Eqns. (12) and noting that terms multiplying the virtual displacements 6v i, Sv 2 
and 64> must vanish identically yields, after eliminating the shear resultant, the equilibrium equations 
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with the associated boundary conditions at s = 0, / 


£vi = 0 

Sv 2=0 (16) 

^ = 0 

where N*, Q*, and M* denote the axial, shear and moment resultants applied at the ends of the arch, 
respectively. 

Expressions for the force and moment resultants, N and M, respectively, are obtained from the 
constitutive law. A unified hereditary visco-elastic-plastic law developed by Walker [22] to characterize 
the time and temperature dependence of Hastelloy X is employed in this study. The selection of Walker’s 
functional theory, a highly generalized representation for a three parameter viscoelastic solid, was based 
principally upon three considerations. First, from a strictly practical point of view, a substantial body of 
experimental work [22-24] had been performed to establish the temperature dependence of the constitu- 
tive law parameters for a wide range of temperatures. These efforts included a validation which examined 
the predictive capability of the law through a comparison of analytical and experimental results for time 
variable thermo-mechanical load cycling of uniaxial specimens. 

A second factor which favored selection of the Walker law was that the law is able to reproduce 
certain forms of classical behavior as a limiting case. For example, saturation of the drag stress parameter 
produces an effect equivalent to isotropic hardening of a material. Similarly, saturation of the backstress 
produces an effect equivalent to kinematic hardening. Finally, the associated laws which govern the 
evolution of the state variables provide for the opportunity to include effects related to both dynamic 
and static thermal recovery. 

The final reason for selection of the Walker law is that it can be expressed in both differential and 
integral formats. For this particular study, the differential format of the law was found to be the most 
convenient. However, it was hope that the availability of an integral format would provide, at a future 
date, the opportunity to extend some prior work involving kinematic bounding of nonlinear integral 
formulations [25] to also include some form of constitutive bounding. 

The general integral form of Walker’s functional theory has been provided in Appendix A. That 
appendix also contains a derivation of the differential form from the integral format. It should be noted 
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that, in the modeling of the elevated temperature behavior of Hastelloy X, additional simplifications in 
the form of the law were possible . These simplifications resulted from the fact that, for Hastelloy X, a 
number of material constants are zero over the entire temperature range. 

As a result of these various simplifications, the differential format of Walker’s functional theory, 
lor one-dimensional loading of Hastelloy X, has the general form 
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Finally, since the time rate of change of temperature is relatively low in the sample problems, 
ail terms where 0 appeared were assumed to be insignificant and therefore neglected. Note that this 
provides some minor simplifications to Eqn. (17b) and makes the reference backstress, OJj, independent 
of the time rate of change in temperature. 

The law is postulated upon an additive decompostion of the total strain into elastic, thermal and 
inelastic components, e et and e c , respectively. Therefore 
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Thus, for the state of one-dimensional loading considered in this study, this yields 
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Therefore, integrating this expression over the cross section and employing the force resultant defini- 
tion, Eqn. (11), yields 
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where E 0 and E x denote the zero and first integrals of the elastic modulus over the cross section 
and Nfl and N c represent “pseudo-resultant” type quantities defined as 
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Note that the first integral of the elastic modulus over the cross section, E\ , does not necessarily vanish 
since it is presumed that the elastic modulus can be a function of temperature which is not necessarily 


constant across a cross section. An interseting aspect of this is that it induces a form of bending- 
stretching coupling {see also Eqn. (22)} similar to that of which can occur in a laminated composite 
material. Finally, it should be noted that, although the quantities defined by Eqn. (21) each have 
the units of a force resultant, these definitions are merely employed as a mathematical convenience to 
simplify subsequent expressions. 

Multiplying Eqn. (19) by the thickness coordinate, 77 , integrating over the cross section and em- 
ploying the moment resultant definition, Eqn. (12), similarly yields 
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where E 2 denotes the second integral of the elastic modulus over the cross section and the “pseudo- 
moment” resultants are defined by 

Me — j 7)Eeo dA and N c = J dA (21) 

Consequently, substituting Eqns. (20) and (22) into Eqns. (15) yields the general governing relations 
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Note that, in the development of these relationships, it is assumed that the moments of the elastic 
modulus, Eq % Ei and E 2 are constant with respect to the axial coordinate. This implies that these 
quantities are independent of stress or strain and also that the temperature is constant along the length 
direction. 

Equations (22) and (23) represent a pair of interrelated spatially dependent nonlinear differential 
relationships which describe the deformation of the beam or arch. In their present form, they are 
stated in terms of the strain along the centroidal axis and the rotation of the cross section-. These two 
■quantities, in turn, are interrelated through the axial and transverse displacements of the centroidal 
axis. Substitution of the appropriate expressions for the centroidal axis strain and cross section rotation 
ultimately results in a set of equations, one of third order and the other of fourth order, in terms of these 
displacement functions. 

Due to the significant nonlinearity of these equations, a numerical method ofsolution was selected. 
The particular method employed is an adaptation of Newton’s method for the solution of nonlinear 
algebraic equations [26]. The basic approach is to develop an iterative procedure in which a “close” 
trial solution can be directed towards the actual solution. Note that this method neither guarantees 
convergence to a solution nor that any solution found is unique. 

The basic methodology is to expand a nonlinear differential term such as dX m dY n , where dX and 
dY represent differentials of the functions X and Y and m and n represent integer exponents, into the 
sum and products of trial solutions X and Y and corrections of the form AX and AY. Consequently, 
substituting X = X + AX and Y = Y + AY yields, 



( 24 ) 
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where /(X, Y)0(AX, AY] represents a nonlinear function of X and Y of second and higher order terms 
in AX and AY. Provided the trial solution is close to the true solution, these higher order terms should 
be small in comparison to the linear terms and thus may be neglected. Consequently, the nonlinear term 
may be closely approximated by the linear form 

dX m dY n % dX m dY n + mdX m - 1 dY n d(AX) + ndX m dY n ” 1 d(AY) (25) 

By employing the technique described above, the nonlinear differential terms appearing in the 
governing equations are approximated in terms of linear differential equations for the corrections to 
assumed trial functions. These coupled differential equations are then converted to a set of coupled 
algebraic relations through the approximation of the linear derivative terms for the corrections to the 
assumed deflections using central difference formulae. Appendix B provides, for example, the general 
finite difference expressions developed for the initially circular arch. 

A matrix iteration procedure is employed to refine an assumed trial solution. Each sucessive set of 
corrections is used to update the trial solution until convergence is obtained. Tests for such convergence 
included the consideration of the magnitude of each set of corrections as well as the overall accuracy for 
which each of the individual nodal equations was solved. In this regard it is noted that an equivalent 
degree of coupling did not exist between the in-plane and transverse equations of equilibrium. Typically, 
the accuracy of the solution of the transverse equation of equilibrium was strongly dependent upon' 
the transverse deflection but only weakly influenced by the in- plane displacement. Conversely, the in- 
plane equation of equilibrium was strongly dependent upon both the transverse and in-plane deflections. 
Consequently, the rate of convergence of the transverse nodal equations was much more rapid than the 
in-plane ones. ' - 

In contrast to the treatment of the equations of equilibrium as a typical boundary value problem, 
the solution for the changes to the moment and force resultants between successive states of quasi-static 
equilibrium is handled as an initial value problem. As such, a fourth order Runge-Kutta integration 
routine was employed to integrate the constitutive law at a preselected set of points across the cross 
section for each axial node used in the finite difference mesh. In this process, it was assumed that the 
change in actual stress at each of these points could be approximated as a linear function of time over 
a given, reasonably short, time interval. 

Once the “trial stresses” at the end of the time interval had been approximated, a Newton-Cotes 
quadrature formula was used to numerically approximate the force and moment resultants. These results 
were then employed to compute the deflections associated with this new state of quasi-static equilibrium. 
From this deflection solution, a revised actual stress field could be calculated and compared with that 
which had been employed to integrate the constitutive law. If more than nominal differences existed 
between the assumed changes in actual stress and those computed, the linear approximating functions 
were adjusted and the entire process repeated. Otherwise, the results were accepted and the analysis 
proceeded on to the next time increment. 



Generally, it was found that after the first few time steps, this procedure worked very well. Most 
of the computation efFort was expended in the solution of the “boundary value” part of the problem and 
not in the repetitive iteration of the material behavior. Principally, this was due to the fact that the rate 
of change in actual stress was reasonably constant and was therefore easy to estimate between successive 
increments. Where this began to break down, however, was when the behavior of the material began to 
resemble a classical plastic response. Under these conditions, the rate of change in actual stress would 
change rapidly even over very short time intervals. Consequently, accurate forecasting its rate of change 
was more difficult and therefore a greater number of iterations was needed to close the numerical loop. 

NUMERICAL RESULTS 

The sample problems considered are an eccentrically loaded cantilever beam-column aiid a simply 
supported pressure loaded shallow circular arch. In both cases, the behavior of the structural element is 
examined for constant loading with constant temperatures of 400, 600 and 800 ° C and with sinusoidal 
temperature variations about the temperatures of 400 and 600 ° C. Simultaneous variations in loading 
and temperature are not examined. However, the potential influences of time-invariant temperature 
gradients in the depth direction are examined for both the constant and variable temperature cases. 

The beam-column considered is 30.48 cm long having a square cross section of 1.27 cm depth 
and thickness. The eccentric load is applied 0.3048 cm below the centroidal axis yielding an eccentricity 
ratio of 0.01. The direction of the load is assumed to remain constant. Twelve axial nodes are used to 
model the beam-column. Additionally, a five point transverse grid is employed at each axial node to 
approximate variations in stress and strain across the cross section. One transverse grid point is located 
at each extreme surface of the cross section and one is positioned on the centroidal axis. The remaining 
two points are spaced equidistantly between these three points. It should be noted that the results from 
this “twelve axial node five grid” model compare favorably with results obtained using greater numbers 
of axial nodes and transverse grid points. 

The circular arch examined is a 8.59 deg. segment of a circle. Physically, the arc length of the 
arch is 22.86 cm with an initial radius of curvature of 152.4 cm. For this relatively shallow arch the rise 
of the centroidal axis is approximately 0.43 cm. Unlike the beam-column, the arch has a rectangular 
“cross section 0.51 cm in width and 0.38 cm in depth. The arch is divided into fourteen segments for the 
numerical model. Again, a five point transverse grid is established at the location of each axial node. 
Similar to the beam-column, the results obtained with this “fourteen- segment five transverse grid model” 
compare favorably with models employing greater numbers of both. Illustrations of the geometry of the 
beam-column and arch models are provided in Figs. 3 and 4. 

At the outset, it is noted that the particular dimensions selected for these sample problems are not 
based upon the consideration of “typical” structural elements. Instead, the dimensions are specifically 
chosen to accelerate the onset of time-dependent behavior. The purpose of this is to minimize the length 
of the initial “stable” response period thereby reducing the overall magnitude of the computational 
effort. In general, a “realistically” sized structural element would be appreciably stiffer and thus provide 
a much longer period of stable response. However, other than this extension of the “stable” useful life, 
the behavioral characteristics of such “realistic” structural elements would be highly similar to those of 
the example problems. 



Finally, before proceeding with a detailed discussion of the numerical results, some introductory 
remarks on the principal factors which determine the form of the response merit consideration. Ex- 
amination of Eqn. (17a), for instance, reveals that the inelastic strain rate is directly determined by 
the relative magnitudes of the actual stress and the backstress. The inelastic strain rate will change 
whenever the rate of change in actual stress varies significantly from the rate of change of the backstress. 
In some situations, the rates of change of actual stress and backstress tend to equilibrate yielding a 
relatively constant difference between these quantities. With this, the rate of inelastic straining tends to 
decrease to a nearly constant value leading to a response which may be characterized as initial “primary 
creep” transitioning to “secondary” or “steady” creep . 

However, in other situations, the rate of change in actual stress increases much more rapidly than 
that of the backstress. Provided the overall magnitude of the difference between them is not excessive, 
this produces a general response akin to that of accelerating or “terciary” creep. The final possibility 
occurs when the difference between the actual stress and the backstress is very large. When this happens, 
the exponential nature of Eqn. (17a) creates a situation where inelastic strain rate tends to follow stress 
increment directly. Consequently, in the limit, the law exhibits a behavior that may be likened to that 
of incremental plasticity. 

Based on the above, it should be evident that the crucial factors in the analysis are those which 
determine how rapidly the actual stress and the backstress change with time. The most significant 
of these factors were found to be a geometrical effect related to the applied bending moment and the 
growth law for the backstress. Consequently, the overall response of the structural element was found 
to be determined by the relative interaction between these effects. 

The geometric effect occurs in both the beam-column and the arch. However, it is more easily 
visualized with respect to the beam-column geometry and therefore is described in that context. Essen- 
tially, the maximum (in magnitude) stresses in the beam-column are determined by the moment created 
by the end load. Transverse defection of the beam-column, increasing the moment arm of the eccentric 
load, increases the magnitude of these stresses. However, tending to counterbalance this effect is the 
end rotation. Because the line of action of the load remains constant, end rotation reduces the effective 
moment arm created by the eccentricity of the load. Thus, this tends to reduce the magnitude of the 
end moment. 

As discussed in a prior study [25], the relative significance of these two opposing influences is 
related to the amount of eccentricity and the time-dependent characteristics of the material. However, 
for an eccentricity of the order employed in this study, these two opposing effects can approximately 
counterbalance one another only when the end deflection is relatively minor. Consequently, beyond 
this threshold, the transverse deflection effect predominates and the end moment inherently increases. 
Therefore, the magnitudes of the maximum stresses increase with deflection and thus also with time. 

In contrast, for slowly varying temperature changes, allowing the 0 terms of Eqn. (17b) to be 
neglected, the evolution of the backstress is governed by a relationship of the form 

(lu=n 2 i e -((l n (20) 

With reference to Eqn. (17e), it should evident that G must be non-negative. Thus, the magnitude 


of the rate of change of the backstress depends upon the signs and magnitudes of the inelastic strain 
rate and the difference between the current and reference values of backstress. Because the actuai 
stress is only indirectly coupled to the backstress, through the elastic and inelastic deformations and 
the constitutive growth laws, an increase in actual stress does not necessarily generate an equivalent 
increase in backstress. Thus, the magnitudes of the quantities appearing in Eqn. (26) need not change in 
equivalent or proportional amounts. As such, the rate of change of backstress may increase, decrease or 
remain relatively constant, depending upon the current values of the various parameters. Consequently, 
this provides a very wide variety of possible results. 

The Beam-Column at Constant Temperature 

The time-dependent deflection of the beam-column at constant temperature provides the simplest 
demonstration of these effects. Figures 5, 6 and 7 illustrate the relative time-dependent end deflection of 
the beam-column under constant load for constant temperatures of 400, 600 and 800 ° C, respectively. 
The magnitude of the loading is expressed in terms of the ratio of the applied load to the Euler load for the 
perfect configuration, P*, where P* = x 2 EI/4L 2 . Note that the Euler load is a function of temperature 
due to the temperature dependence of the elastic modulus. Also note that the relative transverse 
deflection (i.e.: the vertical axis) represents the ratio of the time-dependent transverse deflection to the 
initial elastic deflection. Consequently, these results indicate the relative increase in deflection produced 
by inelastic straining. 

Except for the lowest load at a temperature of 400 ° C, the 400 and 600 ° C beam-column results 
exhibit a short initial settling period followed by a virtually linear increase in transverse deflection 
with respect to time. Generally, this type of behavior is synonomous with a response of primary creep 
transitioning to secondary creep. Not unexpectedly, the higher loadings produce .the greater rates of 
increase. This infers that under these conditions, the difference between the actual stress and the 
backstress must remain nearly constant with the greater numerical differences occurring at the highest 
loadings. 

This hypothesis is, of course, confirmed by the results presented in Fig. 8, which illustrates the 
difference between the maximum (in magnitude) actual stress and the backstress for the 400 °C case.f 
Note that, due to the combination of bending and axial compressive loading, the maximum (in magni- 
tude) actual stress occurs in the extreme fibers adjacent to the wall on the same side of the cenfroida! 
plane as the applied end load. 

Thus since, except for some slight initial variations, the difference between the actual stress and 
backstress remains virtually constant, the right-hand side of Eqn. (17a) effectively is a constant. This 
yields a constant, and in this particular case, a relatively low rate of inelastic straining. Consequently, 
this low rate of inelastic straining does not significantly alter the deflection of the beam-column. As such, 
significant changes in the actual stresses do not occur. Concurrently, the low rate of inelastic straining 

f Results for the difference between the maximum (in magnitude) actual stress and the backstress 
at a 600 ° C constant temperature are not included because they are virually identical to those of Fig. 8 
and are available elswhere [27]. Again, the difference tends to remain relatively constant following a 
short initial period; the higher differences occurring for the higher end loads. 


yields a low the rate of change in bacicstress. Thus, the combined impact of these effects is to maintain 
an approximately constant difference between actual stress and backstress. 

In contrast to the 400 and 600 °C results, virtually all levels of loading at an 800 °C constant 
temperature produce an accelerating rate of transverse deflection. Only the lowest level of loading pro- 
duces a steady creep response; all the other higher loadings result in an accelerating rate of deformation. 
At the two highest load levels, the beam-column deflects so rapidly that it could be considered to have 
failed almost instantaneously. Of course, in comparing results between these different temperature cases, 
the significant differences in the time scales of figures should be kept in mind. At the lower tempera- 
tures, elapsed time can be expressed in hours whereas at the highest temperature, elapsed time must be 
indicated in seconds. 

The reason why the 800 ° C results are so different from those at the lower temperatures is 
directly attributable to the behavior of the difference between the actual stress and the backstress at 
this temperature, illustrated in Fig. 9. Note that, at the lowest load, the initial relaxation is followed by 
a period of constant difference between the actual stress and backstress. Thus, for only this load, the 
increase in actual stress resulting from the increasing transverse deflection effectively is counterbalanced 
by the concurrent growth in the backstress. Higher levels of loading alter the relative rates of growth of 
these quantities. At higher loads, the rate of increase in actual stress far exceeds that of the backstress. 
Consequently, this results in increasing the inelastic strain rate, increasing the rate of deflection and 
thereby further increasing the rate of change of actual stress. The overall process thus reinforces itself, 
thereby accelerating the approach to failure. 

Before proceeding to a discussion of the arch at constant temperature, a few words concerning 
the initial relaxation are warranted. The “relaxation” process is actually a combination of effects. First, 
the inelastic deformation, tends to limit the rate of increase in the actual -stresses at the extreme fibers. 
Sirrudtaneously, due to the requirements for equilibrium, this transfers load bearing responsibility toward 
the centroidal axis. Generally, due to the relatively low initial magnitudes of these stresses, this normally 
does not produce any significant inelastic staining near the centroidal axis. 

The exception to this occurs when the behavior of the beam-column is of the type demonstrated in 
the 800 ° C case. The rapid inelastic straining at the extreme fibers generates significant increases in the 
magnitudes of the stresses throughout the central core of the beam-column. Thus, appreciable inelastic 
straining occurs over the. majority of the cross section. Therefore, in addition to the interaction between 
the maximum (in magnitude) actual stress and the backstress, inelastic straining along the centroidal 
axis begins to strongly influence the overall response. In fact, the form of “failure” which results might 
be characterized as. a viscoelastic analog to a “plastic hinge.” _ 

The Arch at Constant Temperature 

The behavior of the pressure loaded shallow arch at constant temperatures of 400, 600 and 800 ° G 
is generally similar to that of the beam-column. Figures 10 and 11, for example, illustrate the transverse 
deflection at the center of the arch for pressure loadings at the 600 and 800 ° C temperatures. Employing 
the solution method described herein, the critical load is estimated to lie between 179 to 186 kPa (26 to 
27 psi). 

Again, the increase in temperature from 600 to 800 ° C reduces useful life by more than an order 



of magnitude. This is, of course, readily apparent from the significant difference in time scale between 
these two figures. The differences between the actual stress and the backstress for these arch examples 
are similar to the those provided for the beam-column and thus have not been included. 

One common factor to both temperatures is, however, the sensitivity of the results to load mag- 
nitude. Note that only a slight increase in pressure can substantially alter the character of the response. 
The reason for this is that in the arch, unlike the beam-column, the centroidal axis stress is always 
significant and therefore it always influences behavior. Essentially, due to the curvature, the assumed 
boundary conditions and the requirements of equilibrium, significant compressive stessing occurs at the 
centroidal axis. This induces a high rate of compressive inelastic straining along the centroidal axis 
thereby “reducing” the nominal arc length of the arch. This geometric change tends to accomodate 
additional transverse deflection through reduction of both the “nominal” curvature and the “unloaded” 
arc length of the arch. Note also that, unlike the beam-column where such an effect is localized near the 
wall, this centroidal axis straining occurs over a significant part of the central section of the arch. Thus, 
the “failure” zone tends to be distributed and not localized. This is why the sudden rapid increase in 
transverse deflection which occurs with the beam-column is not as pronounced in the arch. 

The Influence of Time-Invariant Temperature Gradients 

The influence of a time-invariant, linear (in the depth direction) temperature gradient on the 
constant temperature response is the next factor considered. This type of temperature gradient is 
typically associated with a steady state constant flow of heat through the depth.* For simplicity, the 
gradient is assumed to remain constant with respect to the length of the element. 

For bookkeeping purposes, temperature gradients for which the temperature is greatest at the 
upper surface and lowest at the bottom surface are considered to be “positive.” Conversely, a situation 
where the greatest temperature is at the bottom surface is denoted as a' “negative” gradient. Finally, 
note that all the gradients were established so that the temperature of the centroidal axis would remain 
at the nominal case temperature of 400, 600 or 800 ° C. Thus, a H-10 9 C gradient for a 400 ® C beam- 
column implies that the upper, centroidal and lower surface temperatures are 405, 400 and 395 ° C, 
respectively. - 

The existence of this type of temperature gradient introduces two effects. First, due to the 
temperature dependence of the elastic modulus, the first integral of the elastic modulus over a cross 
section does not vanish. Thus, this introduces a weak level of bending-stretching coupling into the 
governing equations. The second effect is that such gradients cause thermal bending of the element. 
Positive gradients cause downward bending thereby augmenting the initial mechanically induced bending. 
Conversely, a negative gradient reduces the load induced bending of the element. 

In general, it was found that neither of these effects exert a major influence on the constant 
temperature behavior. Figure 12, for example, illustrates the typical relative increase in transverse 
deflection resulting from a depth direction temperature gradients for the 400 ° C beam-column at a load 
ratio of 0.50. Note that temperature difference a of +10 or —10 °C through the 1.27 cm thick beam- 
column represents a gradient of +7.87 °C/cm. The largest temperature difference considered, namely 
± 25 ° C, corresponds to a heat flux through the depth on the order of 400 kW/m 2 ° C. 

As mentioned above, Fig. 12 provides the relative increase in deflection. In this particular context, 



the relative deflection denotes the ratio of the time-dependent deflection for a given temperature differ- 
ence to the initial elastic deflection of the uniform temperature beam-column. Thus, this ratio indicates 
the gross amount of amplification resulting from inelatic and thermo-elastic efTects. 

In general, temperature differences of 5 ° C and less has little overall impact on the deflection of the 
constant temperature beam-columns. The +10 ° C and larger temperature differences did, however, have 
some impact on the transverse deflection. Typically, gradients of these magnitudes generated observable 
increases in the deflection after several hours of loading. Figure 13 provides a typical example of this 
employing results obtained for a 600 °C beam-column with an applied load ratio of 0.50. The solid line 
represents the initial ratio of end deflection with a thermal gradient to that of the zero gradient case. 
The dashed lined provides the same ratio after 12.5 hours of loading. Note that, in general, although the 
deflection of the beam-column is modified due to the presence of the thermal gradient, the magnitude of 
this change is not sufficient to alter the occurance of the general behaviors illustrated in Figs. 8 through 
11. These results are typical of the behavior observed for the other constant temperature beam-column 
and arch examples which were examined. 

As mentioned above, the imposition of such thermal gradients also creates some weak bending- 
stretching coupling in the governing equations. This also is found to be an insignificant effect. Specifi- 
cally, a number of temperature gradient test cases where the coupling term, E\ t was artifically set to zero 
were examined. There was no appreciable difference in the results obtained from these cases compared 
to results obtained when the coupling term was retained. Consequently, it can be concluded that the 
potential bending-stretching coupling which might be induced by the temperature dependence of the 
elastic modulus is negligible. 

Arch and Beam-Column Variable Temperature Behavior 

The final aspect of behavior examined is the response of the arch and beam-column under si- 
nusoidally varying temperatures. Results are presented for both the arch and beam-column without 
a time-invariant temperature gradient in the depth direction. The combined influences of sinusoidal 
temperature variations and a time-invariant depth direction thermal gradient are investigated only for 
the case of the beam-column. 

In the cases examined the temperature is assumed to^vary in a sinusoidal manner about an elevated 
mean temperature of either 400 or 600 ° C. Amplitudes of 50, 100 and 150 °C are employed in this 
study. Also, in keeping with the utilization of a quasi-static analysis, 1200 and 1800 sec periods are 
employed for the sinusoids. Sinusoidal variations about the 800 °C temperature were not studied due 
to the extremely short useful life which was exhibited in the constant temperature cases. 

Finally, the thermal model includes a slight initial delay between when the load is first applied 
and when the sinusoidal temperature variations commence. Typically, with some combinations of higher 
temperatures and higher loads, a step-like initial transient is generated in the constitutive law state 
variables. Superimposing the start of a sinusoidal temperature at the same time that this step-like 
increase in the state variables occurs* tends to force the Runge-Kutta integration routine to employ 
exceptionally small time step increments to provide accurate results. However, the need to employ 
such small time step increments, which substantially increases the overall computational burden, passes 
almost immediately once the transient change has taken place. Thus, the inclusion of this short delay 


allows the integration routine to stabilize before the start of the temperature oscillations. Typically, a 
delay period on the order of 50 sec is employed. Based on a limited series of test cases, it was found 
that such a short delay period has no appreciable impact on the overall results. 

Figures 14 and 15 illustrate the time dependent deflection of the arch and the beam-column, 
respectively, for 50 and 100 °C amplitude temperature oscillations about a 600 °C temperature. Note 
that the loading employed in .both of these cases would nominally produce a “stable” time-dependent 
response under constant temperature conditions. Several interesting features can be discerned on these 
figures. First, the shape of each deflection curve is a distorted sinusoid. The upper peaks tend to be 
exaggerated whereas the lower peaks are well rounded. Additionally, a distinct underlying increasing 
trend in time-dependent deflection can be observed in both. Specifically, both the upper peak and lower 
peak deflections increase between each cycle. 

The distortion of the peaks from that of a true sinusoid is a result of the temperature dependence 
of the inelastic strain. During the higher temperature portion of the cycle, inelastic straining increases 
due to the temperature sensitivity of the constitutive law parameters. This results in a rapid increase 
in the deflection. Conversely, at lower temperatures, the inelastic strain rate is reduced. Thus, little 
inelastic deformation occurs during this part of the cycle. Consequently, the behavior during the lower 
portion of the temperature cycle closely approximates that of a simple thermo-elastic deformation, where 
deflection would follow the thermal cycle exactly. 

Tables I and II summarize these increasing trends for the arch and beam-column cases, respectively. 
Note that these tables present results obtained for a variety of temperature differences in the depth 
direction. 


Table I. Rate of Change -of Arch Deflection for Sinusoidal Temperature Variations 


Pressure 

(kPa) 

-Sine 

Amplitude 

°C 

Rate of Change 
1 st to 2 nd 
Upper Peak 

(cm/sec) 
1'* to 2 nd 
Lower Peak 

165 

25 

3.18xl0“ 7 

2.74xl0“ 7 

165 

50 

3.84xl0- 7 

3.20xl0“ 7 

165 

100 

4.98x10“° 

2.87x10-° 

152 

100 

7.87xl0“ 7 

6.12xl0- 7 


As Table I illustrates, the magnitude of the increase in deflection between the first and second 
upper peaks appreciably exceeds that between the first and second lower peaks. Although not indicated 
in this or the following table, the increase in deflection at the midpoint temperature of the cycle does 


fall between the values of the extremes. It should be noted that this unequal expansion of the thermal 
cycle deflection indicates that a form of creep ratchetting behavior is occurring. 

This table also demonstrates the significant nonlinearity involved in this phenomena. Note that 
the first doubling of the amplitude of the thermal cycle, from 25 to 50 ° C, only results in a modest 
increase in the peak to peak deflections. However, a second doubling, from 50 to 100 ° C, increases the 
peak to peak deflection by slightly more than an order of magnitude. Table II illustrates that a similar 
effect occurs for the beam-column. Note that except for the 50 ° C amplitude without a temperature 
gradient, increasing the amplitude of the thermal cycle increases the change in peak to peak deflection. 
Again, the increases occur in a nonlinear manner. 

Table II also demonstrates that the increase in deflection between consecutive peaks lias an overall 
decreasing trend. Namely, the absolute increase in deflection decreases as the number of cycles increases. 
The cause of this decreasing trend is attributed to the growth in backstress. Generally, the maximum 
actual stresses do not change significantly between consecutive peaks of the thermal cycle. This is because 
they are principally determined by the thermo-elastic deformation of the element. The accumulated 
inelastic deformation does not significantly alter this stress distribution. However, the amount of inelastic 
strain is sufficient to create an increase in the backstress between consecutive peaks, which exceeds that 
of the actual stress. Consequently, for any particular temperature of the thermal cycle, the difference 
between the actual stress and the backstress decreases. Thus, the rate of inelastic straining is reduced. 
Of course, this effect would ultimately cease to be significant once the backstress approaches saturation. 


Table II. Rate of Change of Beam-Column Deflection for Sinusoidal Temperature Variations 


Sine 

Temp. 


Rate of Change (cm/sec) 


-Ampl. 

Diff. 

Upper Peaks 

Lower Peaks 


°c 

°C 

•r* to 2 nd 

2 nd to y d 

1“ to 2 nd 

2 nd to r d 

50 

0 

1.48x10"* 

1.27x10-® 

1.48x10-* 

1.27x10-* 

50 

• +10 

' 8.46x10-® 

8.26x10-® 

8.05x10-® 

7.82x10"® 

50 

+20 

4.29xl0 -7 

4.06xl0~ 7 

4.01xl0~ 7 

3.84xl0~ 7 

100 

0 • 

1.07x10"® 

1.03x10-® 

9.83xl0" 7 

9.45xl0~ 7 

100 

+ 10 

2.21x10-* 

2.06x10-® 

2.01x10“* 

1.88x10-® 

100 

+20 

4.72x10-® 

4.27x10"® 

4.24x10'® 

3.86x10-® 

150 

- 0 

1.87xl0- 5 

1.66xl0 -5 ‘ 

1.61xl0~* 

1.45x10“® 

150 

+10 

3.22xl0" 5 

2.91xl0- 5 - 

2.77x10-® 

2.52x10"* 

150 

+ 20 

6.41xl0 -5 

4. 93x10“ 5 

4.67x10-* 

4.32x10-* 


The remaining influence illustrated by the data of Table II is the effect that the imposition of 
a time-invariant depth direction temperature gradient has on the thermal cycle inelastic response. In 
general, the change in peak to peak deflection is significantly greater in the presence of the time-invariant 
temperature gradient. Typically, almost an order of magnitude increase in peak to peak change in 


deflection occurs for the lowest temperature cycle amplitude. Interestingly, the effect of the depth 
direction gradient thereafter becomes less significant as the amplitude of the thermal cycle increases. 

The ultimate significance of this effect lies in whether or not temperature gradients in the depth 
direction are truly negligible. The data of this table tend to indicate that the depth direction gradient 
may not be a negligible factor, when the magnitude of the gradient is of the same relative magnitude as 
the thermal cycle. In such cases, the thermal gradient does appear to exert an appreciable influence on 
the overall response. Obviously, this also tends to imply that the impact of the depth direction gradient 
is potentially very significant in transient cases. In such situations, the magnitude of the temperature 
difference through the depth can easily approach that of the gross transient change in temperature. 

CONCLUSIONS 

The elevated temperature behavior of generic types of structural elements fabricated from a typ- 
ical aerospace alloy have been studied analytically by using a nonlinear kinematic analysis to express 
the geometry of deformation and by employing a recently developed nonlinear unified hereditary consti- 
tutive law to express the time-temperature dependence for the material. The results of this study have 
demonstrated that, due to the specific format of the constitutive law, the behavioral response of the 
structural element is determined principally by the difference between the actual stress in the structural 
element and the backstress state variable of the constitutive law. The first of these, the actual stress 
in the element, is basically controlled by the geometry of deformation. However, the second factor, the 
backstress, is governed by the appropriate growth rules of the constitutive law itself. 

This implies that accurate results for such an analysis can be obtained only if the form of the 
backstress growth law and the numerical values employed therein have themselves been established to a 
reasonably high degree of certainty. Otherwise, the prediction of the condifions unde* which the response 
of the structural element may change from that of a “steady stable” form of creep to a rapid approach 
to failure will not be able to be established with any degree of reliability. For entirely similar reasons, 
the actual stresses in the structural element also need to be established accurately to ensure that reliable 
results ensue. 

The specific results obtained from the sample problems indicate that a constitutive law of the type 
proposed by Walker has the capability to model various forms of creep behavior. Under various constant 
temperature and constant load conditions, it has been shown that the predicted response may vary from 
that of simple primary creep-secondary creep to that of almost instantaneous occurring terciary creep. 
Additionally, under varying temperature conditions, the Walker law has demonstrated the capability to 
predict the development of the type of inelastic strain biasing which can produce creep ratchetting. 

Assuming that the Walker law and associated constants provide an accurate representation of 
the elevated temperature response of a material such as Hastelloy X, the sample problem results also 
indicate that the existence of a depth-direction temperature gradient in a thin section may not always 
be negligible. Specifically, under constant temperature conditions, the results indicate that such temper- 
ature gradients do not appear to exert any significant influence on behavioral response. However, such 
gradients do have a significant impact when the overall temperature of the element is not constant. The 
existence of this type of gradient can appreciably accelerate the overall rate of creep. It is indicated that 



this specific effect is most pronounced when the magnitude of the temperature difference through the 
depth of the element is on the order of the amplitude of the gross variation in temperature of the element. 
Consequently, in addition to when such conditions may exist in a “quasi-steady state” situation, this 
implies that such gradients may exert considerable influence in transient thermal problems. 
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Appendix A Constitutive Law 


In its most general integral form, Walker's functional theory has the format: 
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Square bracket terms, of the form /x[0(t)], are used to denote the dependence of the material constants A, 
/x, n°, n, m, n l} n 2 , n 3 , n^, n 5 , n 6 , n 7) K 1} and K 2 on the temperature, 0. However, for clarity, the 
indication of the explicit dependence on temperature will be suppressed in the following development. 

The differential format for Walker’s functional theory is obtained through differentiating the above 
relations with respect to time. Employing Leibniz's rule, differential formats for Eqns. (A-ld), (A-li), 
(A-lg) and (A-ih) are readily obtained. A somewhat simpliried differential format for Eqn. (A-lb) 
may be obtained in the following manner. Differentiating Eqn. (A-lb) with respect to time and using 
Leibniz’s rule, after rearranging, yields 


= rif. + (., + n,)o u +0,^0 + [^6 - »,G] jT ' ^ L (.4 - 2) 


The integral which appears in the above is identical to the integral of Eqn. (A-lb). Consequently, 
solving Eqn. (A-lb) for the integral and substituting into Eqn. (A-2) yields the differential format given 
by Eqn. (A-6b). 

With some manipulation, the differential form of Eqn. (A-la) may be stated in a “power law” type 
format. First, Eqn. (A-la) is differentiated with respect to time and then simplified using the deviatoric 
stress tensor to yield 

hh- ««> - i f 2 ^ - P» TT - f ^ i «• (- 1 - 3 > 

Note that to establish the above it is necessary to show that H** = = 0. These, however, follow 

directly from Eqns. (A-lb) and (A-le) provided the inelastic portion of the deformation may be treated 
as, at least, approximately incompressible (i.e., c** « 0). 

Therefore, substituting Eqn. (A-3) into Eqn. (A-2) and using the differential form of Eqn. (A-ld) 

yields: 
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The Q term can be replaced with the aid of the differential forms of Eqns. (A-lf) and (A-lh). After 
some algebraic manipulation and minor rearranging this yields: 
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Consequently, substituting Eqn. (A-5) back into the differential form of Eqn. (A-lf) and using Eqn. (A-4) 
yields, after some algebra, the differential power law of Eqn. (A-la). Thus, the complete set of differential 
relations for Walker’s functional theory are: 
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Note that a differential form for the function Q is not required since it does not appear in any of the 
other expressions. 

Numerical constants for the law were established [22-24] from uniaxial bar specimens tested 
under fully reversed, strain-rate controlled cyclic stress-strain tests. These tests were conducted for 
a variety of temperatures and strain rates. The testing conditions were sufficiently rigorous to cause 
plastic deformation during the loading cycle. The creep and relaxation properties were deduced from 
observation of the behavior of the samples when the cyclic loading was “held” at various points on the 
“steady-state” hysteresis loop. Appendix Figs. A-l through A-8 illustrate the temperature dependent 
parametric values established under that work. Note that the material constants K 2 , ni, n 4 , ns, and n 
are zero over the enitre temperature range considered. This results in additional simplifications to the 
constitutive law. 


Appendix B Finite Difference Equations 


Solutions for the governing nonlinear differential equations are obtained using a Newton type 
method for nonlinear differential equations. It is assumed that trial solutions for the transverse and 
axial deflection functions, denoted as w and u t respectively, exist. From these, trial values for cen- 
troidal axis strain and cross section rotation, i 0 and <£, respectively, are determined. Then, using the 
approximations e 0 = e 0 + Ae 0 and (f> = 4> + A <f> to substitute into Eqns. (16) yields 
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The terms N and M are evaluated by substituting e 0 and 4> into Eqns. (11). Numerical methods are used 
to approximate the derivative terms on the right-hand side of Eqns (B-l). The second derivative of the 
moment correction, which appears on the left-hand side of Eqn. (B-la), also is evaluated numerically. 
This avoids the numerical problems which can result in approximating fourth derivatives directly. 

Assuming w = w + Aw and a similar expansion for u, upon substitution into Eqn. (7) and 
neglecting terms of order (Aw) 2 and higher, yields 
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where the factors A and B are defined as: 
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Differentiating Eqn. (B-2) with respect to s yields, after rearranging, 
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where C and D are given by 
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Developing an expression for A<£ is accomplished by using Eqn. (8). Since it is assumed that <f> = + A </>, 

substitution into Eqn. (8) yields 
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where the approximations sin A (f> « A<£ and cos A</> ss 1 have been used. For small e Q it is reasonable 
to assume that Ae a < e 0 , thus allowing the Ae a term to be neglected. Thus, since the first term on the 
left-hand side, according to Eqn, (8), is equal to the first term of the right-hand side, this implies that 
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Therefore, combining Eqns. (8) and (B-7) yields 
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Derivatives of A(/> are obtained from Eqn. (B-8). In this, it is assumed that the term 1/A may be treated 
as an approximate constant and thus need not be differentiated. This provides significant simplifications 
without engendering any substantial inaccuracies since the correction terms become negligible as the 
exact solution is approached. Consequently, differentiating Eqn. (B-8) with respect to s yields: 
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Central difference formulae are used to approximate the derivative terms. However, prior to 
substituting into Eqns. (B-l) it is useful to establish the following definitions: 
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Thus, substituting the various expressions given above into Eqns. (B-l) yields, after rearranging, 
the general finite difference nodal equations 
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These two equations provide the finite-difference approximation for the deflection of the internal nodes 
of the beam, that U for the subscript i over the range 3 < i < nr— 1, (see Fig. 4). Therefore, for an arch ol 
n + i total nodes, they provide a system of 2n — 6 equations in 2n -f 2 unknowns. Thus, eight additional 
equations are required to provide a unique solution for the problem. Six of these equations are obtained 
by using the boundary condition. These provide two sets of three equations, one set applicable to node l 
and the other set for node n + 1. 

As evidenced above, the boundary conditions do not provide a sufficient number of relations to 
enable unique solution. This, in part, results from the inherent coupling generated by the assumption of 
Euler-Bernoulli bending when both axial and transverse deflections are possible. Note that the deflection 



functions are not only related to one another, but they also appear in the same functional format in 
both expressions. The derivative of one is added to (subtracted from) the other multiplied by the 
curvature. This, in turn, indicates that a form of implicit coupling exists between the two generalized 
displacements, centroidal axis strain and rotation. Ideally, the generalized displacements should be 
completely independent. 

A number of approaches can be employed to generate the two additional equations needed for 
unique solution. The one used in this study is to require the “centroidal axis” strain in the wall to 
vanish. Mathematically, this is equivalent to appending, to the existing system of equations, the two 
additional equations: 


Ae 0 \ - —Z 0 \ and Ae or = -£ or , (B - 14a, fc) 

where e Q \ and e OT represent the strain in the left and right extensions, respectively. No condition is 
established for the “rotations” which might occur in the wall sections. This is not felt to be significant 
since the nodal mesh, at least for the equilibrium problem, is a centroidal mesh. 

Using the boundary conditions and the two “wall strain” relations it is possible to eliminate the 
transverse and axial displacement components for nodes 1, 2, n, and n + 1 from the general finite- 
difference equations. This leaves a system of exactly 2n — 6 equations in the same number of unknowns. 
If the coefficients of these equations are arranged in matrix form, the result is a banded matrix with a 
bandwidth of six. 
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Figure 1. Geometry of Deformation. 
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2 External Forces and Moments Acting on the Arch. 






Figure 6 Beam-Column Deflection at 600 °C. 



Figure 7 Beam-Column Deflection at 800 °C. 





0 400 800 

TIME, sec 


Figure 9 Difference Between Actual Stress and Backstress at 8( 
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Figure 10 Center Deflection Versus Time for the 600 °C Low Arch. 



Figure 11 Center Deflection Versus Time for the 800 3 C Low Arch. 
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Figure 14 Lotv Arch Deflection for 50 and 100 "C Temperature Variations 


P / P g * 0. 5 0 


200 "C pk-pk 


100 *C pk-pk 


TIME, sec 


Figure 15 End Deflection for 50 and 100 °C Amplitude Temperature Variations 












Figure A-3 n -1 as a Function of Temperature 



T 7 ; as a Function of Temperature 


Figure A-4 
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Figure A- 7 Reference Bacxstress. Cl c 33 a Function of Temperature 



